	subroutine wave_prop
	include 'cst.inc'
c****************************************************************************************************
c    This subroutine is related to the wave biased bedload-dominanted transport
c    In Niedoroda et al. (1995), the wave propagation in 1-D is perpandicular to the shore
c    In the 2-D case here, the wave propoagation direction is allowed to have a angle with 
c    the on-shore direction, and the angle is aw 
c 
c    determines the local 'time-averaged' wave propagation direction
c    (aw) and the corresponding breaker angle (ab) by approximating
c    wave refraction and shoaling. the scheme assumes
c    that wave refraction will cause the breaker angle to be reduced
c    from the angle of the offshore waves (acent) with the shoreline (sangle).
c    the amount of this reduction depends 'in some way' on the relative
c    steepness of the offshore shallow bathymetry.  for the present 
c    purposes we assume that the reduction of this angle is constant over
c    the region and can simply be controlled by a factor (rfact) which will
c    have a values between 0 and 1. this is read in as an input. 
c
c    it is further assumed that refraction does not make much difference in
c    water depths greater than 30m. between this depth and the time-averaged
c    breaker depth (taken to be 3m) refraction acts as an exponential 
c    function of water depth (see the function below). these assumptions
c    determine the wave propagation direction (aw) for each grid cell.
c
c    the effect of refraction in changing the wave energy density because of
c    converging or diverging wave rays is also approximated. for present 
c    purposes, this is assumed to vary linearly with the second derivative
c    of the shoreline curve (map view). a rough scaling argument suggests that
c    for most cases edfact is about 0.5. this factor, like rfact, will vary
c    with the steepnest of the offshore profile in some, yet to be determined
c    manner.
c************************************************************************************************
c
c --- Set the breaker angle ---
c     The scheme assumes that wave refraction will cause the breaker angle (ab) to be reduced
c     from the angle of the offshore waves (acent) with the shoreline (sangle). for the present 
c     purposes we assume that the reduction of this angle is constant over the region and can 
c     simply be controlled by a factor (rfact) which will have a values between 0 and 1. 
c     this is read in as an input.
c
c     The description below is related to the concept of power equivalent wave 
c	Wave has directions. For all measured wave directions, they may fall into a cloud. 
c     The principal direction of cloud is called acent. It is the angle of wave crest with
c     a straight (east-west) shoreline. When the shoreline is not straignt, acent is adjusted
c     by the shoreline angle (sangle)      
c
c     For each shoreline element, there is a breaker angle
	do i=1,imax
      	 ab(i)=rfact*(acent-sangle(i)) 
      enddo	 
c
c --- Set the wave propagation direction for each cell ---
c     it is further assumed that refraction does not make much difference in
c     water depths greater than 30m. between this depth and the time-averaged
c     breaker depth (taken to be 3m) refraction acts as an exponential 
c     function of water depth (see the function below). these assumptions
c     determine the wave propagation direction (aw) for each grid cell.
c
	do i=1,imax
	   jj=jshore(i)
	   do j=1,jj
	      if(h(i,j).lt.hwave_max)then !water depth greater than 30m, determined using equation of wave length in page 163 of Komar 
	        aw(i,j)=acent
	      else
	        if(h(i,j).gt.hwave_min)then
	          aw(i,j)=ab(i)
	        else                   !The exponential function is entirely empirical
        	    aw(i,j)=acent-(exp(0.15*(h(i,j)-have_min))*(acent-ab(i)))
              endif  
	      endif    
         enddo
      enddo
c 	write(*,*)'w_prpr  ab'
c	write(*,*)(ab(i),i=1,imax)
c	pause


c	write(*,*)'w_prpr  aw'
c	write(*,*)(aw(3,j),j=1,jshore(3))
c	pause
c
c
c    set the factor that approximates the change of wave energy density at
c    the time-averaged breaker line (shoal) all along the shore. edfact
c    is a factor that is input at the beginning. 
c 
c    the effect of refraction in changing the wave energy density because of
c    converging or diverging wave rays is also approximated. for present 
c    purposes, this is assumed to vary linearly with the second derivative
c    of the shoreline curve (map view). a rough scaling argument suggests that
c    for most cases edfact is about 0.5. this factor, like rfact, will vary
c    with the steepnest of the offshore profile in some, yet to be determined
c    manner.
	ii=imax-1
	do i=2,ii
	   dydx1=(yshore(i)-yshore(i-1))/(xshore(i)-xshore(i-1))
	   dydx2=(yshore(i+1)-yshore(i))/(xshore(i+1)-xshore(i))	
	   dy2dx2=(dydx1-dydx2)/((xshore(i+1)-xshore(i-1))/2)
   	   shoal(i)=(1.0+(edfact*dy2dx2))
      enddo
	shoal(1)=shoal(2)
	shoal(imax)=shoal(ii)

	return
	end